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Summary. The seminal work by Bardeen, Cooper and Schrieffer (BCS) extended 
further by Eliashberg to the intermediate coupling regime solved one of the major 
scientific problems of Condensed Matter Physics in the last century. The BCS theory 
provides qualitative and in many cases quantitative descriptions of low-temperature 
superconducting metals and their alloys, and some novel high-temperature super- 
conductors like magnesium diboride. The theory has been extended by us to the 
strong-coupling regime where carriers are small lattice polarons and bipolarons. 
Here I review the multi-polaron strong-coupling theory of superconductivity. At- 
tractive electron correlations, prerequisite to any superconductivity, are caused by 
an almost unretarded electron-phonon (e-ph) interaction sufficient to overcome the 
direct Coulomb repulsion in this regime. Low energy physics is that of small po- 
larons and bipolarons, which are real-space electron (hole) pairs dressed by phonons. 
They are itinerant quasiparticles existing in the Bloch states at temperatures below 
the characteristic phonon frequency. Since there is almost no retardation (i.e. no 
Tolmachev-Morel- Anderson logarithm) reducing the Coulomb repulsion, e-ph inter- 
actions should be relatively strong to overcome the direct Coulomb repulsion, so 
carriers must be polaronic to form pairs in novel superconductors. I identify the 
long-range Frohlich electron-phonon interaction as the most essential for pairing in 
superconducting cuprates. A number of key observations have been predicted or 
explained with polarons and bipolarons including unusual isotope effects and upper 
critical fields, normal state (pseudo)gaps and kinetic properties, normal state dia- 
magnetism, and giant proximity effects. These and many other observations provide 
strong evidence for a novel state of electronic matter in layered cuprates, which is a 
charged Bose-liquid of small mobile bipolarons. 



1 Introduction 

While a single polaron problem has been actively researched for a long time 
(for reviews see ref. Q] II II El H H IS ED] and the present volume.), multi- 
polaron physics has gained particular attention in the last two decades. For 
weak electron-phonon coupling, A < 1, and the adiabatic limit, uj/Ep <C 1), 
Migdal theory describes electron dynamics in the normal Fermi-liquid state 
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[TTj . and BCS-Eliashberg theory in the superconducting state [TH Q3] (here 
and further I use % = c = ks = 1)- While the electron-phonon (e-ph) interac- 
tion is weak Migdal's theorem is perfectly applied. The theorem proves that 
the contribution of diagrams with "crossing" phonon lines (so called "ver- 
tex" corrections) is small if the parameter Xw/Ep is small, where A is the 
dimensionless (BCS) e-ph coupling constant, u> is the characteristic phonon 
frequency, and Ef is the Fermi energy. Neglecting the vertex corrections, 
Migdal [11] calculated the renormalized electron mass as m* = mo(l + A) 
(near the Fermi level), where mo is the band mass in the absence of e-ph in- 
teraction, and Eliashberg [13] extended Migdal's theory to describe the BCS 
superconducting state at intermediate values of A < 1. Later on many re- 
searchers applied Migdal-Eliashberg theory with A even larger than 1 (see, 
for example, Ref. [14], and references therein). With increasing strength of 
interaction and increasing phonon frequency, w, finite bandwidth [TB] and 
vertex corrections |17] become increasingly important. But unexpectedly for 
many researchers who applied the non-crossing approximation even at A > 1 
we have found that the Migdal-BCS-Eliashberg theory (with or without ver- 
tex corrections) breaks down entirely at A ~ 1 for any value of the adiabatic 
ratio lo/Ef since the bandwidth is narrowed and the Fermi energy, Ep is 
renormalised down exponentially so the effective parameter Xco/Ep becomes 
large [18]. 

The electron-phonon coupling constant A is about the ratio of the electron- 
phonon interaction energy E p to the half bandwidth D w N(Ep)' 1 , where 
N(E) is the density of electron states in a rigid lattice. One expects [18] that 
when the coupling is strong, A > 1, all electrons in the bare Bloch band are 
"dressed" by phonons since their kinetic energy (< D) is small compared 
with the potential energy due to the local lattice deformation, E p , caused by 
electrons. In this strong coupling regime the canonical Lang-Firsov transfor- 
mation pjj] can be used to determine the properties of the system. Under cer- 
tain conditions [501 [HI HH H3] , the multi-polaron system is metallic but with 
polaronic carriers rather than bare electrons. This regime is beyond Migdal- 
Eliashberg theory, where the effective mass approximation is used and the 
electron bandwidth is infinite. In particular, the small polaron regime cannot 
be reached by summation of the standard Feynman-Dyson perturbation dia- 
grams using a translation-invariant Green's function G(r, r', r) = G(r — r', r) 
with the Fourier transform G(k, fl) prior to solving the Dyson equations on 
a discrete lattice. This assumption excludes the possibility of local violation 
of the translational symmetry [24] due to the lattice deformation in any or- 
der of the Feynman-Dyson perturbation theory similar to the absence of the 
anomalous (Bogoliubov) averages in any order of perturbation theory [11] . To 
enable electrons to relax into the lowest polaronic band, one has to introduce 
an infinitesimal translation-noninvariant potential, which should be set zero 
only in the final solution obtained by the summation of Feynman diagrams for 
the Fourier transform G(k, k',.!?) of G(r, r',r) rather than for G(k, Q) [15j . 
As in the case of the off-diagonal superconducting order parameter, the off- 
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diagonal terms of the Green function, in particular the Umklapp terms with 
k' = k + G, drive the system into a small polaron ground state at sufficiently 
large coupling. Setting the translation-noninvariant potential to zero in the 
solution of the equations of motion restores the translation symmetry but in 
a polaron band rather than in the bare electron band, which turns out to be 
an excited state. 

Alternatively, one can work with momentum eigenstates throughout the 
whole coupling region, but taking into account the finite electron bandwidth 
from the very beginning. In recent years many numerical, and analytical stud- 
ies have confirmed the conclusion [18] , that the Migdal-Eliashberg theory 
breaks down at A > 1 (see, for example Refs. [25l [26l H3 IM [Ml ffl [311 
|32l [331 O [35l [361 EH [38] and contributions to this book). With increasing 
phonon frequency the range of validity of the 1/ A polaron expansion extends 
to smaller values of A [35]. As a result, the region of applicability of the Migdal- 
Eliashberg approach (even with vertex corrections) shrinks to smaller values 
of the coupling, A < 1, with increasing to. Strong correlations between carriers 
reduce this region further (see [30]). 

Carriers in the fascinating advanced materials are strongly coupled with 
high-frequency optical phonons, making small polarons and non-adiabatic ef- 
fects relevant for high-temperature superconductivity, colossal magnetoresis- 
tance phenomenon, and molecular electronic devices (see Part III). Indeed the 
characteristic phonon energies 0.05 — 0.2 eV in cuprates, manganites and in 
many organic materials are of the same order as generally accepted values of 
the hopping integrals 0.1 — 0.3 eV. 

As reviewed in this book lattice polarons are different from ordinary elec- 
trons in many aspects, but perhaps one of the most remarkable difference is 
found in their superconducting properties. By extending the BCS theory to- 
wards the strong interaction between electrons and ion vibrations, a charged 
Bose gas (CBG) of tightly bound small bipolarons was predicted by us [41] 
with a further prediction that high critical temperature T c is found in the 
crossover region of the e-ph interaction strength from the BCS-like polaronic 
to bipolaronic superconductivity [18j . This contribution describes what hap- 
pens to the conventional BCS theory when the electron-phonon coupling be- 
comes strong. The author's particular view of cuprates is also presented. 

2 Electron-phonon and Coulomb interactions in Wannier 
representation 

For doped semiconductors and metals with a strong electron-phonon interac- 
tion it is convenient to transform the Bloch states |k) to the site (Wannier) 
states |m) using the canonical linear transformation of the electron operators, 




(1) 
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where i = (m, s) includes both site m and spin s quantum numbers, and N is 
the number of sites in a crystal. In the site representation the electron kinetic 
energy takes the following form 

H e = [T(m - m')<W - fiSij] c\c , (2) 

where 

k 

is the bare hopping integral in the rigid lattice, /i is the chemical potential, 
j = (n, s'), and is the bare Bloch band dispersion. 

The electron-phonon and Coulomb interactions acquire simple forms in 
the Wannier representation, if their matrix elements in the momentum repre- 
sentation 7(q,f) and V c (q) depend only on the momentum transfer q, 



H e -ph = ^i^i Mq,")<V + H - c -] ' ( 3 ) 

and 



H e -e = ^ ^ y c (m - n)hihj. (4) 



2 . 

1+3 



Here 



and 



uM = _J= 7 ( q ,„) e *i-> (5) 

Vc(m) = ±J2 V ^y tl ' m > ( fi ) 
q 

are the matrix elements of e-ph and Coulomb interactions, respectively, in 
the Wannier representation for electrons, hi = c\ci is the electron density 
operator, and <i q „ annihilates the ^-branch phonon with the wave vector q 
and frequency w q „. Taking the interaction matrix elements depending only 
on the momentum transfer one neglects terms in the electron-phonon and 
Coulomb interactions, which are proportional to the overlap integrals of the 
Wannier orbitals on different sites. This approximation is justified for narrow 
band materials with the bandwidth 2D less than the characteristic value of 
the crystal field. As a result, the generic Hamiltonian takes the following form 
in the Wannier representation, 

H = Y \- T ( m ~ m ')tiss' - fJ^ij] 4 C 3 + X] w q""» t u »(<l> I/ ) < V + H.c] 

~ V c (m - n)h t h 3 + J2 ^ v {d\ v d^ v + 1/2). (7) 



2 . 
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Here we confine our discussions to a single electron band and the e-ph matrix 
element depending only on the momentum transfer q. This approximation 
allows for qualitative and in many cases quantitative descriptions of essen- 
tial polaronic effects in advanced materials. There are might be degenerate 
atomic orbitals in solids coupled to local molecular-type Jahn- Teller distor- 
tions, where one has to consider multi-band electron energy structures (see 
[42]). 

The quantitative calculation of the matrix element in the whole region of 
momenta has to be performed from pseodopotentials [43, 44, 45 . On the other 
hand one can parameterize the e-ph interaction rather than to compute it from 
first principles in many physically important cases [46j . There are three most 
important interactions in doped semiconductors, which are polar coupling 
to optical phonons (i.e the Frohlich e-ph interaction), deformation potential 
coupling to acoustical phonons, and the local (Holstein) e-ph interaction with 
molecular type vibrations in complex lattices. While the matrix element is 
ill defined in metals, it is well defined in doped semiconductors, which have 
their parent dielectric compounds, together with bare phonons w q „ and the 
electron band structure Here the effect of carriers on the crystal field and 
on the dynamic matrix is small while the carrier density is much less than 
the atomic one (for phonon self-energies and frequency renormalizations in 
polaronic systems see Ref. [20l [47] ) . Hence one can use the band structure 
and the crystal field of parent insulators to calculate the matrix element in 
doped semiconductors. 

The e-ph matrix element 7(q) has different g-dependence for different 
phonon branches. In the long wavelength limit (q <C n/a, a is the lattice 
constant), 7(q) cx q n , where n = —1, and n = — 1/2 for polar optical, molec- 
ular (u>q = ujq)) and acoustic ((w q cx q)) phonons, respectively. Not only q 
dependence is known but also the absolute values of 7(q) are well parameter- 
ized in this limit. For example in polar semiconductors the interaction of two 
slow electrons at some distance r is found as (see below) 

"W = ^w-^Ew , i)iV qT (8) 

q 

The Coulomb repulsion in a rigid lattice is V c (r) = e 2 /e oc r, and the second 
term represents the difference between the Coulomb repulsion screened with 
the core electrons and the repulsion screened with both core electrons and 
ions. Hence the matrix element of the Frohlich interaction depends only on 
the dielectric constants and the optical phonon frequency u)$ as 

l7(q)| 2 = — , 0) 

KU>o 

where K = (e" 1 - e^ 1 ) -1 - 

One can transform the e-ph interaction further using the site-representation 
also for phonons. The site representation of phonons is particularly convenient 
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for the interaction with dispersionless local modes, whose w„ v = u) v and the 
polarization vectors e qi/ = e„ are q independent. Introducing the phonon 
site-operators 

rf ™ = 4=5> <q ' n eV ( 10 ) 
V-iV k 

one transforms the deformation energy and the e-ph interaction as [23j 

H ph = ^u v {dl v d Ilv + l/2), (11) 

and 

+ dn„), (12) 

n,m,i/ 

respectively. 

Here g u {m) is a dimensionless force acting between the electron on site 
m and the displacement of ion n, and e m _ n = (m — n)/|m — n| is the unit 
vector in the direction from the electron m to the ion n. This real space 
representation is convenient in modelling the clcctron-phonon interaction in 
complex lattices. Atomic orbitals of an ion adiabatically follow its motion. 
Therefore the electron does not interact with the displacement of the ion, 
whose orbital it occupies, that is g v (0) = 0. 



3 Breakdown of Migdal-Eliashberg theory in the 
strong-coupling regime 

Obviously a perturbative approach to the e-ph interaction fails when A > 1. 
However one might expect that the self-consistent Migdal-Eliashberg (ME) 
theory is still valid in the strong-coupling regime because it sums the infinite 
set of particular non-crossing diagrams in the electron self-energy. One of the 
problems with such an extension of the ME theory is a lattice instability. The 
same theory applied to phonons yields the renormalised phonon frequency Q — 
2A) 1 / 2 [H]. The frequency turns out to be zero at A = 0.5. Because of this 
lattice instability Migdal [11] and Eliashberg [13] restricted the applicability 
of their approach to A < 1. However, it was shown later that there was no 
lattice instability, but only a small renormalisation of the phonon frequencies 
of the order of the adiabatic ratio, u)/fj, <C 1, for any value of A, if the adiabatic 
Born-Oppenheimer approach was properly applied [48j . The conclusion was 
that the Frohlich Hamiltonian correctly describes the electron self-energy for 
any value of A, but it should not be applied to further renormalise phonons. 

In fact, ME theory cannot be applied at A > 1 for the reason, which 
has nothing to do with the lattice instability. Actually the 1/A muti-polaron 
expansion technique [20] shows that the many-electron system collapses into 
the small polaron (or bipolaron) regime at Aw 1 for any adiabatic ratio. 
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To illustrate the point let us compare the Migdal solution of the simple 
molecular-chain Holstein model 49] with the exact solution [55] in the adi- 
abatic limit, u>o/t — > , where t = T(a) is the nearest neighbour hopping 
integral. The Hamiltonian of the model is 



Ei 1 d 2 kx 2 
' ^2Mdx 2 + ~2~ 



(13) 



where Xi is the normal coordinate of the molecule (site) i, and k = Mw . The 
Migdal theorem is exact in this limit. Hence in the framework of the Migdal- 
Eliashberg theory one would expect the Fermi-liquid behaviour above T c and 
the BCS ground state below T c at any value of A. In fact, the exact ground 
state is a self-trapped insulator at any filling of the band, if A > 1. 

First we consider a two-site case (zero dimensional limit), i,j — 1,2 with 
one electron, and than generalise the result for an infinite lattice with many 
electrons. The transformation X = (xi+xz), £ = X\ — xi allows us to eliminate 
the coordinate X , which is coupled only with the total density [n\ +ri2 = 1). 
That leaves the following Hamiltonian to be solved in the extreme adiabatic 
limit M — > oo: 



H 



-t(c\c 2 + 4ci) + {\kt) 1/2 i(c\ Cl - c\c 2 ) 



ke 



(14) 



The solution is 
where 

a 

and the energy is 



ij = (ac\+(3cl)\0), 
t 



[t 2 + ({Xkt) 1 / 2 ^ + (t 2 + Afct^ 2 ) 1 /2)2]l/2 ' 

(Afct^g+^ + Afc^ 2 ) 172 
[t 2 + ((\kty/ 2 Z + it 2 + Afci£ 2 ) 1/2 ) 2 ] 1/2 ' 



E=^~(t 2 + \kte) 1/2 - 



(15) 
(16) 
(17) 

(18) 



In the extreme adiabatic limit the displacement £ is classical, so the ground 
state energy Eq and the ground state displacement £o are obtained by min- 
imising Eq.(18) with respect to £. If A > 0.5 we obtain 



and 



'f(4A 2 - 1) 
Xk 



1/2 



(19) 
(20) 
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Fig. 1. The ground-state energy (in units of t, solid line) and the order parameter 
(thin solid line) of the adiabatic Holstein model. The Migdal solution is shown as 
the dashed line 



The symmetry-breaking "order" parameter is 

^ = ^-^= [2A+(4A2 - 1)1/2]2 ~ 1 (21) 
1 [2A + (4A2-1)V2]2 + 1 - V 1 ) 

If A < 0.5, the ground state is translation invariant, A = 0, and -Bo = —t, £ = 
0, ft = —a. Precisely This state is the "Migdal" solution of the Holstein model, 
which is symmetric (translation invariant) with \a\ = \(3\. When A < 0.5, the 
Migdal solution is the only solution. However, when A > 0.5 this solution is 
not the ground state of the system, Fig.l. The system collapses into a localised 
adiabatic polaron trapped on the "right" or on the "left-hand" site due to the 
finite lattice deformation £o ^ 0. 

The generalisation for a multi-polaron system on the infinite lattice of 
any dimension is straightforward in the extreme adiabatic regime. The adi- 
abatic solution of the infinite one-dimensional (ID) chain with one electron 
was obtained by Rashba [50j in a continuous (i.e. effective mass) approxima- 
tion, and by Holstein [49] and Kabanov and Mashtakov [25[ [47] for a discrete 
lattice. The latter authors studied the Holstein two-dimensional (2D) and 
three-dimensional (3D) lattices in the adiabatic limit. According to Ref. [25] 
the self-trapping of a single electron occurs at A > 0.875 and at A > 0.92 in 
2D and 3D, respectively. The radius of the self-trapped adiabatic polaron, r p , 
is readily derived from its continuous wave function [50j 

t/>(r) ~ l/cosh(Ar/a). (22) 

It becomes smaller than the lattice constant, r p = a/X for A > 1. Hence 
the multi-polaron system remains in the self-trapped insulating state in the 
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strong-coupling adiabatic regime, no matter how many polarons it has. The 
only instability which might occur in this regime is the formation of self- 
trapped bipolarons, if the on-site attractive interaction, 2\zt, is larger than 
the repulsive Hubbard U [FT]. Self-trapped on-site bipolarons form a charge 
ordered state due to a weak repulsion between them [41j[36] (see also [52]). 

The transition into the self-trapped state due to a broken translational 
symmetry is expected at 0.5 < A < 1.3 (depending on the lattice dimen- 
sionality) for any electron-phonon interaction conserving the on-site electron 
occupation numbers. For example, Hiramoto and Toyozawa 153"] calculated the 
strength of the deformation potential, which transforms electrons into small 
polarons and bipolarons. They found that the transition of two electrons into 
a self-trapped small bipolaron occurs at the electron-acoustic phonon coupling 
A ~ 0.5, that is half of the critical value of A at which the transition of the 
electron into the small acoustic polaron takes place in the extreme adiabatic 
limit, lu << zt (here z is the coordination lattice number). The effect of the 
adiabatic ratio cu/zt on the critical value of A was found to be negligible. The 
radius of the acoustic polaron and bipolaron is about the lattice constant, 
so that the critical value of A does not very much depend on the number of 
electrons in this case either. As discussed below the non-adiabatic corrections 
(phonons) allow polarons and bipolarons to propagate as the Bloch states in 
narrow bands. 



4 Polaron dynamics 



4.1 Polaron band 



A self-consistent approach to the multy-polaron problem is possible with the 
"1/ A" expansion technique [2U] , which treats the kinetic energy as a perturba- 
tion. The technique is based on the fact, known for a long time, that there is 
an analytical exact solution of a single polaron problem in the strong-coupling 
limit A — * oo [19 . Following Lang and Firsov we apply the canonical trans- 
formation e s to diagonalise the Hamiltonian. The diagonalisation is exact, if 
T(m) = (or A = oo): 

H = e s He~ s , (23) 



where 



S = - 22 hi [ui(q,v)dq V - H.c] 



(24) 



is such that — —S. The electron and phonon operators are transformed as 



q = Ci exp 



y^ttf(q,i/)d qi/ - H.c. 



(25) 



and 
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<V ^ d ^-Yl ™i u i( ( l' L ')' ( 26 ) 

i 

respectively. It follows from Eq.(27) that the Lang-Firsov canonical transfor- 
mation shifts the ions to new equilibrium positions. In a more general sense 
it changes the boson vacuum. As a result, the transformed Hamiltonian takes 
the following form 



& = Eft* - ^M c j - E pJ2 hi+ Yl ^vK^ + V 2 ) + |l> 

(27) 



2 / j -ij n i n ji 

j i n,v ijij 



where 



Gij = T(m - n)6 ss > exp ^^[tij(q,i>) - u»(q,i')]<V - H.c^j (28) 

is the renormalised hopping integral depending on the phonon operators, and 
Vij = v(m — n) = (29) 
V c (m l7(q^)| 2 ^ cos [q ' ( m - n )] 

is the interaction of polarons comprising their Coulomb repulsion and the 
interaction via a local lattice deformation. In the extreme infinite-coupling 
limit, A — > oo, we can neglect the hopping term of the transformed Hamilto- 
nian. The rest has analytically determined eigenstates and eigenvalues. The 
eigenstates \N) = |nj,n q „) are sorted by the polaron n ms and phonon n q „ 
occupation numbers. The energy levels are 

E = -( / u + £ l p )y^n» + ^ ^ VijUiUj + ^ w q „(ra q „ + 1/2), (30) 

% ijij q 

where rij = 0, 1 and n q „ = 0, 1, 2, 3, ....oo. 

The Hamiltonian, Eq.(27), in zero order with respect to the hopping de- 
scribes localised polarons and independent phonons, which are vibrations of 
ions relative to new equilibrium positions depending on the polaron occupa- 
tion numbers. The phonon frequencies remain unchanged in this limit. The 
middle of the electron band falls by the polaron level shift E p due to a poten- 
tial well created by lattice deformation, Fig. 2, 

^^EhMlV ( 31 ) 

Now let us discuss the 1/A expansion. First we restrict the discussion to 
a single-polaron problem with no polaron-polaron interaction and /i = 0. The 
finite hopping term leads to the polaron tunnelling because of degeneracy of 
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the zero order Hamiltonian with respect to the site position of the polaron. 
To see how the tunnelling occurs we apply the perturbation theory using 1/ A 
as a small parameter, where 

A=§, (32) 

and D = zt. The proper Bloch set of TV-degenerate zero order eigenstates with 
the lowest energy {—E p ) of the unperturbed Hamiltonian is 

i k '°> = 7^E c - cx p( ik - m )i°>< ( 33 ) 

* m 

where |0) is the vacuum. By applying the textbook perturbation theory one 
readily calculates the perturbed energy levels. Up to the second order in the 
hopping integral they are given by 

E(k) = -E p + e k (34) 

KME^Alk> q ,)| 2 

where |k', n q „) are the exited states of the unperturbed Hamiltonian with one 
electron and at least one real phonon. The second term in Eq.(34), which 
is linear with respect to the bare hopping T(m), describes the polaron-band 
dispersion, 

e k = ^T(m)e-s 2 ( m )cxpHk-m), (35) 



where 



9 2 



(m) = -^-$>(q^)| 2 [l-cos(q.m)] ( 36 ) 



2N 



is the band-narrowing factor at zero temperature. The third term in Eq.(34), 
quadratic in T(m) , yields a negative almost k- independent correction to the 
polaron level shift of the order of 1/A 2 . The origin of this correction, which 
could be much larger than the fist-order contribution (Eq.(35) contains a small 
exponent), is understood in Fig. 2. The polaron localised in the potential well 
of the depth E p on the site m, hops onto a neighbouring site n with no 
deformation around and comes back. As any second order correction this 
transition shifts the energy down by an amount of about —t 2 /E p . It has little 
to do with the polaron effective mass and the polaron tunneling mobility 
because the lattice deformation around m does not follow the electron. The 
electron hops back and forth many times (about e 9 ) waiting for a sufficient 
lattice deformation to appear around the site n. Only after the deformation 
around n is created does the polaron tunnel onto the next site together with 
the deformation. 

Oddly enough, analysing the Holstcin two-site model some authors took 
the second-order correction in t in Eq. (34) as a measure of the polaron motion 
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E=0 



Fig. 2. "Back and forth" virtual transitions of the polaron without any transfer of 
the lattice deformation from one site to another. These transitions shift the middle 
of the band further down without any real charge derealization 



and arrived at an erroneous conclusion that "polarons are no longer describ- 
able in terms of quasiparticles having a well-defined dispersion" [54) and "the 
Lang-Firsov approach, which is generally believed to become exact in the limit 
of antiadiabaticity and an electron-phonon coupling going to infinity, actually 
diverges most from the exact results precisely in this limit..." [55, 56J. Later on 
it became clear [ST] [55] 023 that this controversy is the result of the erroneous 
identification of the polaron kinetic energy by the authors of Refs. [54l l55l [56] . 

4.2 Polaron spectral and Green's functions 

The multi-polaron problem has an exact solution in the extreme infinite- 
coupling limit, A = oo, for any type of e-ph interaction conserving the on- 
site occupation numbers of electrons. For the finite coupling 1/A perturbation 
expansion is applied. The expansion parameter is actually [HI ED HO] 

1 

« 1, 



2zA 2 



so that the analytical perturbation theory has a wider region of applicabil- 
ity than one can expect using a semiclassical estimate E p > D. However, 
the expansion convergency is different for different e-ph interactions. Exact 
numerical diagonalisations of vibrating clusters, variational calculations (see 
Refs. [2BH17lliai271|31|37l|35]and[B3|53]), dynamical mean-field approach 
in infinite dimensions [34], and Quantum-Monte-Carlo simulations (see Refs. 
[62[3i[65l[6a[6Zl[6a[6i[2n] and [H [71]) simulations revealed that the 



Superconducting Polarons and Bipolarons 



13 



ground state energy (w — E p ) is not very sensitive to the parameters. On 
the contrary, the effective mass, the bandwidth, and the polaron density of 
states strongly depend on the adiabatic ratio ioft and on the radius of the 
interaction. The first-order in 1/A perturbation theory is practically exact in 
the non-adiabatic regime lo > t for any value of the coupling constant and 
any type of e-ph interaction. However, it overestimates the polaron mass by 
a few orders of magnitude in the adiabatic case, uj <C t, if the interaction is 
short-ranged [26] . 

A much lower effective mass of the adiabatic Holstcin polaron compared 
with that estimated using the first order perturbation theory is the result 
of poor convergency of the perturbation expansion owing to the double-well 
potential [49 in the adiabatic limit. The tunnelling probability is extremely 
sensitive to the shape of this potential and also to the phonon frequency dis- 
persion. The latter leads to a much lighter Holstein polaron compared with 
the nondispersive approximation [72] . Importantly, the analytical perturba- 
tion theory becomes practically exact in a wider range of the adiabatic pa- 
rameter and of the coupling constant for the long-range Frohlich interaction 
[38]. 

Keeping this in mind, let us calculate the one-particle GF in the first 
order in 1/A. Applying the canonical transformation we write the transformed 
Hamiltonian as 

H = Hp + H p h + Hi nt , (37) 

where 

iJ p = ^(k) C t Ck (38) 

k 

is the "free" polaron contribution, 

H ph = £w q (dtd q + l/2) (39) 
q 

is the free phonon part (spin and phonon branch quantum numbers are 
dropped here), and £k = Z'E^ — u is the renormalised polaron-band dis- 
persion. The chemical potential fi includes the polaron level shift — E p , and 
it could also include all higher orders in 1/ A corrections to the polaron spec- 
trum, which are independent of k. The band-narrowing factor Z' is defined 
as 

7 , = E m T(m) e -g 2 ^exp(-zk-m) 

£ m T(m)expHk.m) ' ^ 

which is Z' — exp(—-fEp/uj)) with 7 < 1 depending on the range of the e- 
ph interaction and phonon frequency dispersions. The interaction term Hi nt 
comprises the polaron-polaron interaction, Eq.(29), and the residual polaron- 
phonon interaction 

H P- P h = ~ &j) ph ]4 c 3, ( 41 ) 
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where (<Jij) ph means averaging with respect to the bare phonon distribution. 
We can neglect H p _ p f l in the first-order of 1/A <C 1. To understand spectral 
properties of a single polaron we also neglect the polaron-polaron interaction. 
Then the energy levels are 

Efh = Z* n * + + V2], (42) 



and the transformed eigenstates \fh) are sorted by the polaron Bloch-state 
occupation numbers, nk = 0, 1, and the phonon occupation numbers, n q = 
0,1,2, ...,oo. 

The spectral function of any system described by quantum numbers m, n 
and eigenvalues E n , E m is defined as (see, for example [25] ) 

A(k, u>) = tt(1 + e-"/ T )e n ' T £ e~ E ^ T \ (n\ c k |m) | 2 6(tu nm + u). (43) 

It is real and positive, A(k, u>) > 0, and obeys the important sum rule 

1 f 00 

- / dujA(k,uj) = 1. (44) 

7T ./-oo 

Here i? is the thermodynamic potential. The matrix elements of the electron 
operators can be written as 



^^e-^nlc^lm) (45) 



[n\ c k \m) = _ 

by the use of the Wannier representation and the Lang-Firsov transformation. 
Here 



Xi = exp 



^u 4 (q)d q - H.c. 



Now, applying the Fourier transform of the (^-function in Eq.(43), 

5(u nm +u) = — dte^-+^, 
2tt J-oo 

the spectral function is expressed as 
1 f°° 1 

A(k, W ) = - / dte^- ^- m) x (46) 

~°° in . n 

{((c i (t)l i (t)ctxj)} + ((c J tlt c . W l. (t) ^}. 

Here the quantum and statistical averages are performed for independent 
polarons and phonons, therefore 
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[ciWMQxM)) = ((a(t)c])) ({Mt)x})) . (47) 

The Heisenberg free-polaron operator evolves with time as 

<*(*) = c k e-^ 4 , (48) 



and 



= jf E ^ (k '- m - k "- n) ((c„(*)4,)) = (49) 

k',k" 

± £[1 - n(k')]e ik '-< m - n >-** t ) 

k' 

c ^(*))) = l^^kOe^''^-")-^'* (50) 



TV 

k' 

where n(k) = [1 + exp^k/^ 1 ] -1 is the Fermi-Dirac distribution function of 
polarons. The Heisenberg free-phonon operator evolves in a similar way, 

dq(t) = <i q e- iw «*, 

and 

= J] ((exp[ Ui (q,t)rf q - tf.c.] exp[- Uj (q)d q - H.c])) , (51) 

q 

where Uij(q,t) = u it} ;(q)e _lWc| *. This average is calculated using the operator 
identity 

e A + 6 = e A e & e-^/\ (52) 

which is applied for any two operators A and B, whose commutator [A, B] 
is a number. Because [rf q ,d q ] = 1, we can apply this identity in Eq.(51) to 
obtain 

e [« i (q,t)d q -Jf.c.] e [-u 3 (q)d q -i?.c] _ e (a*d q -ad q ) x 

e [« i ( q ,t)«*( q )-«*( q ,t)u j ( q )]/2 

where a = Wj(q,i) — Uj(q). Applying once again the same identity yields 

g[u i (q,i)d q -Jf.c.]g[-u.,(q)d q -.Ef.c.] _ £ u* d q(= -Q<i qe -|a| 2 /2 x ^ggj 

e [ui(q,t)uJ(q)-uJ(q,t)«^(q)]/2_ 

Now the quantum and statistical averages are calculated by the use of 

^ e a*d q g-ad q \^ = e -\»\ 2 n^ ^ 

where n u = [exp(aj q /T) — 1] _1 is the Bose-Einstein distribution function of 
phonons. Collecting all multiplies in Eq.(51) we arrive at 
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X < (i)Xj))=e X p|- 5 L^| 7 (q)| 2 / q (nx-n,i)|, 



(55) 



where 



/ q (m,i) = [1 — cos(q • m) cos(w q i)] coth — ^ + i cos(q • m) sin(w q i). (56) 

Here we used the symmetry of "/(— q) = 7(q), because of which terms contain- 
ing sin(q • m) disappeared. The average ^X}Xi(t)^ , which is a multiplier 
in the second term in the brackets of Eq. (46) , is obtained as 



X\Xi(t) 



Xi{t)X) 



(57) 



To proceed with the analytical results we consider low temperatures, T <C w q , 
when coth(aj q /2T) w 1. Then expanding the exponent in Eq.(55) yields 



x i (t)^)) = z y £ 



{E q l7(q) 



2gi[q-(m— n)— o) q t] 



1=0 



(2N) l n 



(58) 



where 



Z = cxp 



(59) 



Then performing summation over m. n, k' and integration over time in Eq.(46) 
we arrive at [75] 



A{k, uj) = j2 [ A t } ( k : «o + 4 +) ( k > W ) 



(=0 



(60) 



where 



A ; (_) (k,w) = ttZ E 



IlLl l7(qr) 
(27V)';! 



(61) 



1 — n j k— q r 



r=l 



and 



4 +) (k, W ) = 7rZ £ 



qi,...qi 



nUl7(qr)| ; 
(2JV)'l! 



n 



r=l 



(62) 
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Clearly Eq.(60) is in the form of a perturbative multi-phonon expansion. Each 
contribution A\ f (k, u>) to the spectral function describes the transition from 
the initial state k of the polaron band to the final state k± Y^ r =i Ir- with the 
emission (or absorption) of I phonons. 

The 1/A expansion result , Eq.(60), is applied to low-energy polaron ex- 
citations in the strong-coupling limit. In the case of the long-range Frohlich 
interaction with high-frequency phonons it is also applied in the intermediate 
regime [38, 69 . Different from the conventional ME spectral function there 
is no imaginary part of the self-energy since the exponentially small at low 
temperatures polaronic damping [23] is neglected. Instead the e-ph coupling 
leads to the coherent dressing of electrons by phonons. The dressing can be 
seen as the phonon "side-bands" with I > 1. 

While the major sum rule, Eq.(44) is satisfied, 



the higher-momentum integrals, J_ duju> p A(]i., to) with p > 0, calculated us- 
ing Eq.(60), differ from the exact values (see Part III) by an amount propor- 
tional to 1/A. The difference is due to a partial "undressing" of high-energy 
excitations in the side-bands, which is beyond the first order 1/A expansion. 

The spectral function of the polaronic carriers comprises two different 
parts. The first (I — 0) k-dependent coherent term arises from the polaron 
band tunnelling, 



The spectral weight of the coherent part is suppressed as Z <C 1. However 
in the case of the Frohlich interaction the effective mass is less enhanced, 
£k = Z'E k - /x, because Z « Z' < 1. 

The second incoherent part ^4i nco / l (k, lu) comprises all the terms with I > 
1. It describes the excitations accompanied by emission and absorption of 
phonons. We notice that its spectral density spreads over a wide energy range 
of about twice the polaron level shift E p , which might be larger than the 
unrenormalised bandwidth 2D in the rigid lattice without phonons. On the 
contrary, the coherent part shows a dispersion only in the energy window of 
the order of the polaron bandwidth, w — Z'D. It is interesting that there is 
some k dependence of the incoherent background as well [73 , if the matrix 
element of the e-ph interaction and/or phonon frequencies depend on q. Only 
in the Holstein model with the short-range dispersionless e-ph interaction 
7(q) =70 and o-> q = luq the incoherent part is momentum independent (see 




(63) 



A coh (k,u)= 4,- J (k, w )+A^(k,a;) 



ttZS(uj - £ k ). 



(64) 



also Ref. 05]), 
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Z °° -y 21 

A incoh (k,u)=w-J2^x (65) 



N ^ 2H 



1=1 



£ {[1 - n (k')] S (w-iwo-f^) + n (k') 5 (w+Jwo-^)} ■ 

k' 

As soon as we know the spectral function, different Green's functions (GF) 
are readily obtained using their analytical properties, indexpolaron! Green's 
function For example, the temperature GF is given by the integral ^23J 

SM-f^. (66) 

where LOk — nT(2k + 1), k = 0, ±1, ±2, .... Calculating the integral with the 
spectral density Eq.(65) we find in the Holstein model [71] 



Here the first term describes the coherent tunnelling in the narrow polaron 
band while the second k-independent sum describes the phonon-side bands. 



5 Polaron-polaron interaction and small bipolaron 

Polarons interact with each other, Eq.(29). The range of the deformation sur- 
rounding the Frbhlich polarons is quite large, and their deformation fields are 
overlapped at finite density. Taking into account both the long-range attrac- 
tion of polarons owing to the lattice deformations and their direct Coulomb 
repulsion, the residual long-range interaction turns out rather weak and re- 
pulsive in ionic crystals pQ. The Fourier component of the polaron-polaron 
interaction, u(q), comprising the direct Coulomb repulsion and the attraction 
mediated by phonons is 

4-7TP 2 

*(q) = — a-Mq)rW (68) 
£oo<r 

In the long- wave limit (q <C 7r/a) the Frohlich interaction [75] dominates in 
the attractive part, which is described by Eq.(9). Fourier transforming Eq.(68) 
yields the repulsive interaction in real space, 

e 2 

u(m-n) = — ; r > 0. (69) 

e |m - n| 

We see that optical phonons nearly nullify the bare Coulomb repulsion in ionic 
solids, where eo S> 1, but cannot overscreen it at large distances. 



Superconducting Polarons and Bipolarons 



19 



Considering the polaron-phonon interaction in the multi-polaron system 
we have to take into account dynamic properties of the polaron response func- 
tion [59j . One may erroneously believe that the long-range Frohlich interaction 
becomes a short-range (Holstein) one due to the screening of ions by heavy 
polaronic carriers. In fact, small polarons cannot screen high-frequency optical 
vibrations because their renormalised plasma frequency is comparable with or 
even less than the phonon frequency. In the absence of bipolarons (see below) 
we can apply the ordinary bubble approximation to calculate the dielectric 
response function of polarons at the frequency fl as 



e(q, Q) 



1 



2,(q)V" (k + q) '" (k) 



(70) 



This expression describes the response of small polarons to any external field 
of the frequency fl < u>o, when phonons in the polaron cloud follow the 
polaron motion. In the static limit we obtain the usual Debye screening at 
large distances (q — > 0). For the temperature larger than the polaron half- 
bandwidth, T > ui, we can approximate the polaron distribution function 
as 

n ( _ (2 - n)e k ' 
2a 3 V 2T 



and obtain 



rc(k) 



(71) 



c(«,0) = l + 2|, 



(72) 



where 



Is = 



27re 2 n(2 - n) 



e Q Ta 3 



1/2 



and n is the number of polarons per unit cell. For a finite but rather low 
frequency, loq > fl 3> w, the polaron response becomes dynamic, 



e(q, Q) - 1 



n 2 



where 



w p(q) = 2«(q) "( k )( £ k+ q - £k) 



k 



(73) 



(74) 



is the temperature-dependent polaron plasma frequency squared. The polaron 
plasma frequency is rather low due to the large static dielectric constant, 
e 3> 1, and the enhanced polaron mass m* 3> m e . 

Now replacing the bare electron-phonon interaction vertex 7(q) by a 
screened one, 7 sc (q, w ), as shown in Fig. 3, we obtain 



7sc(q,Wo) 



7(q) 



s 7(q) 



(75) 
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Fig. 3. E-ph vertex, 7(q) screened by the polaron-polaron interaction, v (q) (dashed- 
dotted line) . Solid and dotted lines are polaron and phonon propagators, respectively. 



because u)q > uj p . Therefore, the singular behaviour of 7(q) ~ 1/q is unaffected 
by screening. Polarons are too slow to screen high-frequency crystal field oscil- 
lations. As a result, the strong interaction with high-frequency optical phonons 
in ionic solids remains unscreened at any density of small polarons. 

Another important point is a possibility of the Wigner crystallization of 
the polaronic liquid. Because the net long-range repulsion is relatively weak, 
a relevant dimensionless parameter r s — m*e 2 /eo(47rn/3) 1 / 3 is not very large 
in ionic semiconductors. The Wigner crystallization appears around r s ~ 100 
or larger, which corresponds to the atomic density of polarons n < 10 -6 with 
£o = 30 and to* = 5m e . This estimate tells us that polaronic carriers are 
usually in the liquid state at relevant doping levels. 

At large distances polarons repel each other. Nevertheless two large po- 
larons can be bound into a large bipolaron by an exchange interaction even 
with no additional e-ph interaction but the Frohlich one [76j [21 1114] . When a 
short-range deformation potential and molecular-type (i.e. Holstein) e-ph in- 
teractions are taken into account together with the Frohlich interaction, they 
overcome the Coulomb repulsion at a short distance of about the lattice con- 
stant. Then, owing to a narrow band, two small polarons easily form a bound 
state, i.e. a small bipolaron. Let us estimate the coupling constant A and 
the adiabatic ratio wo/t, at which the "bipolaronic" instability occurs. The 
characteristic attractive potential is \v\ = D/(X — fi c ), where fi c is the dimen- 
sionless Coulomb repulsion, and A includes the interaction with all phonon 
branches. The radius of the potential is about a. In three dimensions a bound 
state of two attractive particles appears, if 
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\v\> 



Substituting the polaron mass, m* = [2a 2 t] 1 exp(7A-D/<x>o)j we find 



ma* 

24 



U 



4z(A - fi c ) 



(76) 



(77) 



As a result, small bipolarons form at A > [i c + 7r 2 /4z almost independent of 
the adiabatic ratio. In the case of the Frohlich interaction there is no sharp 
transition between small and large polarons, and the first-order 1/A expansion 
is accurate in the whole region of the e-ph coupling, if the adiabatic parameter 
is not too small. Hence we can say that in the antiadiabatic and intermediate 
regime the carriers are small polarons independent of the value of A if the 
e-ph interaction is long-ranged. It means that they tunnel together with the 
entire phonon cloud no matter how "thin" the cloud is. 



6 Polaronic superconductivity 

The polaron-polaron interaction is the sum of two large contributions of the 
opposite sign, Eq.(29). It is generally larger than the polaron bandwidth and 
the polaron Fermi-energy, ep = Z'Ep. This condition is opposite to the weak- 
coupling BCS regime, where the Fermi energy is the largest. However, there is 
still a narrow window of parameters, where bipolarons are extended enough, 
and pairs of two small polarons are overlapped similar to the Cooper pairs. 
Here the BCS approach is applied to nonadiabatic carriers with a nonretarded 
attraction, so that bipolarons are the Cooper pairs formed by two small po- 
larons [18j . The size of the bipolaron is estimated as 

n * ww^> (78) 

where A is the binding energy of the order of an attraction potential v. The 
BCS approach is applied if r^ 3> n -1 / 3 , which puts a severe constraint on the 
value of the attraction 

M « e F . (79) 

There is no "Tolmachev" logarithm in the case of nonadiabatic carriers, be- 
cause the attraction is nonretarded as soon as €f < ^o- Hence a supercon- 
ducting state of small polarons is possible only if A > /i c . This consideration 
leaves a rather narrow crossover region from the normal polaron Fermi liquid 
to a superconductor, where one can still apply the BCS mean-field approach, 

< A - n c -C Z' < 1. (80) 

In the case of the Frohlich interaction Z' is about 0.1-^0.3 for typical values of 
A. Hence the region, Eq.(80), is on the border-line from the polaronic normal 
metal to a bipolaronic superconductor (section 8). 
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In the crossover region polarons behave like fermions in a narrow band 
with a weak nonretarded attraction. As long as A ^> l/\/2z, we can neglect 
their residual interaction with phonons in the transformed Hamiltonian, 



(81) 



written in the Wannier representation . If the condition Eq.(80) is satisfied, 
we can treat the polaron-polaron interaction approximately by the use of the 
BCS theory. For simplicity let us keep only the on-site v and the nearest- 
neighbour v i interactions. At least one of them should be attractive to ensure 
that the ground state is superconducting. Introducing two order parameters 

A) = -v (c m:T c m!i ), (82) 

A\ = --Ul(c m ,TCm+a,l) (83) 

and transforming to the k-space results in the familiar BCS Hamiltonian, 

H p = £k4ck s + ]T[A4 T cL ki + H.c], (84) 

k,s k 

where £ k = £k — M is the renormalised kinetic energy and 

A k = A -A 1 ^±^ (85) 
w 

is the order parameter. 

The Bogoliubov anomalous averages are found as 

(c kT c kX ) = tanh ^ + ^ , (86) 

and two coupled equations for the on-site and inter-site order parameters are 

A) = ~E ,f" tanh ^I±^, (87) 

A 1 = — — V" ^& + "> tanh ^ + 3. (88) 

These equations are equivalent to a single BCS equation for zi k = ^(£k)> 
but with the half polaron bandwidth w cutting the integral, rather than the 
Debyc temperature, 



r^ 11 am 

A(H)= dr,N p ( V )V{S,ri) \" tanh 



2T 



(89) 
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Here Vfarj) = -v - zvi(£ + n)(ri + fi)/w . 

The critical temperature T c of the polaronic superconductor is determined 
by two linearised equations in the limit Aq i — » 0, 

,2 ' 



l + A 



^4 

w 



EL 

w 2 



w 



0. 



(l + B)A 1 = 0, 



(90) 
(91) 



where A = Aq — Ai/s/w, and 



A = 



B = 



ZVi 

2w 

ZV\ 

~2w 



W—fJ, 

— W — fJ, 
W—fl 

-ID — fJ, 



dr) 
dq 



tanh 
V 

r\ tanh ^r- 



These equations are applied only if the polaron-polaron coupling is weak, 
Iuq i| < w. A nontrivial solution is found at 



Lb 

\.\Aw\j 1 exp 



2w 



(92) 



if wo + zvifi 2 /w < 0, so that superconductivity exists even in the case of 
the on-site repulsion, vq > 0, if this repulsion is less than the total intersite 
attraction, z\vi\. There is a nontrivial dependence of T c on doping. With a 
constant density of states in the polaron band, the Fermi energy ep sa ll is 
expressed via the number of polarons per atom n as 



jit = w(n — 1), 



so that 



1.14w-\/ n(2 — n) exp ( - 
V 1 



2w 



(93) 



(94) 



+ zvi[n - l] 2 / 

which has two maxima as a function of n separated by a deep minimum in the 
half-filled band (n = 1), where the nearest-neighbour contributions to pairing 
are canceled. 



7 Mobile small bipolarons 

The attractive energy of two small polarons is generally larger than the po- 
laron bandwidth, A — /i c 3> Z' . When this condition is satisfied, small bipo- 
larons are not overlapped. Consideration of particular lattice structures shows 
that small bipolarons are mobile even when the electron-phonon coupling is 
strong and the bipolaron binding energy is large [401 170| . Here we encounter 
a novel electronic state of matter, a charged Bose liquid, qualitatively differ- 
ent from the normal Fermi-liquid and the BCS superfluid. The Bosc-liquid is 
stable because bipolarons repel each other (see below). 
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7.1 On-site bipolarons and bipolaronic Hamiltonian 

The small parameter, Z'/(X — /i c ) -C 1, allows for a consistent treatment of 
bipolaronic systems [T5], HJJ • Under this condition the hopping term in the 
transformed Hamiltonian H is a small perturbation of the ground state of 
immobile bipolarons and free phonons, 

H = H a + H pert , (95) 

where 

Hq = VijclctcjCi + ^V^cV + 1/2] (96) 

i,j I," 

and 

H pert = ^%c| C j (97) 

Let us first discuss the dynamics of onsite bipolarons, which are the ground 
state of the system with the Holstein non-dispersive e-ph interaction [5TJ [UJ . 
The onsite bipolaron is formed if 

2E p > U, (98) 

where U is the onsite Coulomb correlation energy (the Hubbard U). The 
intersite polaron-polaron interaction is just the Coulomb repulsion since the 
phonon mediated attraction between two polarons on different sites is zero in 
the Holstein model. Two or more onsite bipolarons as well as three or more 
polarons cannot occupy the same site because of the Pauli exclusion principle. 
Hence, bipolarons repel single polarons and each other. Their binding energy, 
A = 2E p — U, is larger than the polaron half-bandwidth, A 3> w, so that 
there are no unbound polarons in the ground state. H pert , Eq.(97), destroys 
bipolarons in the first order. Hence it has no diagonal matrix elements. Then 
the bipolaron dynamics, including superconductivity, is described by the use 
of a new canonical transformation exp(S*2) 41J, which eliminates the first 
order of H pert , 

( S i)fv-l^ E E ■ (99) 

i,j 1 r 

Here E^ p and |/), \p) are the energy levels and the eigenstates of H - Neglect- 
ing the terms of the order higher than (w/A) 2 we obtain 

(H b ) ff ,= (e s 'He- s ') ffi , (100) 

{ H b)ff> » (Ho) f f> (f\°w c \ c i'\p)Wor c \ c Af) x 

1 1 



E p - Ef E p - E f 
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5*2 couples a localised onsite bipolaron and a state of two unbound polarons 
on different sites. The expression (100) determines the matrix elements of 
the transformed bipolaronic Hamiltonian H b in the subspace |/), \ f) with 
no single (unbound) polarons. On the other hand, the intermediate bra (p\ 
and ket \p) refer to configurations involving two unpaired polarons and any 
number of phonons. Hence we have 



where are phonon occupation numbers (0, 1, 2, 3...oo). This equation is an 
explicit definition of the bipolaron binding energy A which takes into account 
the residual inter-site repulsion between bipolarons and between two unpaired 
polarons. The lowest eigenstates of Hb are in the subspace, which has only 
doubly occupied cj ns c| ns , |0) or empty |0) sites. On-site bipolaron tunnelling 
is a two-step transition. It takes place via a single polaron tunneling to a 
neighbouring site. The subsequent tunnelling of its "partner" to the same site 
restores the initial energy state of the system. There are no real phonons 
emitted or absorbed because the (bi)polaron band is narrow. Hence we can 
average Hb with respect to phonons. Replacing the energy denominators in 
the second term in Eq.(lOO) by the integrals with respect to time, 




(101) 




we obtain 




(102) 




Here $™ m / (t) is a multiphonon correlator. 




(103) 



Xj(t) and Xi'(t) commute for any j(q_,v) = j(— <l,v). xj and Xj> commute 
as well, so that we can write 



xj(t)X t ,(t) = J^ e [«i'(q,t)-Mq,t)]*.--ff.c.] ) ( 10 4) 



q 



xjXj, = JJ e [»j'W- u i( < 0]«»Q-ff-c.] ) ( 105 ) 



q 



where the phonon branch index v is dropped for transparency. Applying twice 
the identity Eq.(52) yields 
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x}{t)X v {t)x]x r = IJe^ie-^e-l"! 2 / 2 



(106) 



e [« i /(q,t)-u i (q,t)][u*,(q)-u*(q)]/2-H.c. 

where 

/3 = Ui(q,t) - Mj'(q,t) + Uj(q) - Uj(q). 
Finally using the average Eq.(54) we find 



where 



= p -a 2 (m-m') p -g 2 (n-n') x 



,(t) = e 



(107) 



^(m, m', n, n') = cos[q • (n' — m)] + cos[q • (n — m')] — (108) 
cos[q • (n' — m')] — cos[q • (n — m)]. 

Taking into account that there are only bipolarons in the subspace, where 
Hb operates, we finally rewrite the Hamiltonian in terms of the creation 6^ = 
c mt c m| an( i annihilation b m = c m |C m | operators of singlet pairs as 



-E 

m 

E 



t(m — m'jbjjjbm' + ^w(m — m')r 



(109) 



There are no triplet pairs in the Holstein model, because the Pauli exclusion 
principle does not allow two electrons with the same spin to occupy the same 
site. Here n m = bl n b m is the bipolaron site-occupation operator, 



v(m — m') = 4u(m — m') + tr ' (m — m'), 



(110) 



is the bipolaron-bipolaron interaction including the direct polaron-polaron 
interaction t>(m — m') and a second order in T(m) repulsive correlation 



v (2 >(m-m ; ) = 2i / dte~ iAt <P™ r 
Jo 



m 
mm 



(t). 



(in) 



This additional repulsion appears because a virtual hop of one of two polarons 
of the pair is forbidden, if the neighbouring site is occupied by another pair. 
The bipolaron transfer integral is of the second order in T(m) 



t(m-m') = -2iT 2 (m-m') / dte~ lAt $™™,(t). 
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The bipolaronic Hamiltonian, Eq.(109) describes the low-energy physics of 
strongly coupled electrons and phonons. Using the explicit form of the multi- 
phonon correlator, Eq.(107), we obtain for dispersionless phonons at T <C w , 



(S mm 
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Expanding the time-dependent exponents in the Fourier series and calculating 
the integrals in Eqs.(112) and (111) yield [78] 
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When A -C Wo, we can keep the first term only with Z = in the bipolaron 
hopping integral, Eq.(113). In this case the bipolaron half-bandwidth zt(&) is 
of the order of 2w 2 /(zA). However, if the bipolaron binding energy is large, 
A ^> wo, the bipolaron bandwidth dramatically decreases proportionally to 
e~ 4g in the limit A — > oo. However, this limit is not realistic because A = 
2E p — V c < 2g 2 ujQ. In a more realistic regime, ui < A < 2g 2 ujQ, Eq.(113) 
yields 
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On the contrary, the bipolaron-bipolaron repulsion, Eq.(114) has no small ex- 
ponent in the limit A — * oo, oc D 2 /A. Together with the direct Coulomb 
repulsion the second order i/ 2 ) ensures stability of the bipolaronic liquid 
against clustering. 

The high temperature behavior of the bipolaron bandwidth is just the 
opposite to that of the small polaron bandwidth. While the polaron band 
collapses with increasing temperature [S], the bipolaron band becomes wider 

EH], 



i(m) oc —= exp 



E„ + A 



2T 



(116) 



for T > ujq. 



7.2 Superlight intersite bipolarons in the Frohlich-Coulomb model 
(FCM) 

Any realistic theory of doped ionic insulators must include both the long- 
range Coulomb repulsion between carriers and the strong long-range electron- 
phonon interaction. From a theoretical standpoint, the inclusion of the long- 
range Coulomb repulsion is critical in ensuring that the carriers would not 
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form clusters. Indeed, in order to form stable bipolarons, the e-ph interaction 
has to be strong enough to overcome the Coulomb repulsion at short dis- 
tances. Since the realistic e-ph interaction is long-ranged, there is a potential 
possibility for clustering. The inclusion of the Coulomb repulsion V c makes 
the clusters unstable. More precisely, there is a certain window of V c /E p in- 
side which the clusters are unstable but mobile bipolarons form nonetheless. 
In this parameter window bipolarons repel each other and propagate in a 
narrow band. 

Let us consider a generic "Frohlich-Coulomb" Hamiltonian, which explic- 
itly includes the infinite-range Coulomb and electron-phonon interactions, in 
a particular lattice structure [80]. The implicitly present infinite Hubbard U 
prohibits double occupancy and removes the need to distinguish the fermionic 
spin, as soon as we are interested in the charge excitations alone. Introducing 
spinless fermion operators c n and phonon operators d mu , the Hamiltonian is 
written as 

H = J2 T ( n - n ')4cn' + \ E - n')cUnclc n , + (117) 




The e-ph term is written in the real space representation (section 2), which is 
more convenient in working with complex lattices. 

In general, the many-body model Eq.(117) is of considerable complexity. 
However, we are interested in the non/near adiabatic limit of the strong e-ph 
interaction. In this case, the kinetic energy is a perturbation and the model 
can be grossly simplified using the Lang-Firsov canonical transformation in 
the Wannier representation for electrons and phonons, 

S = 2J .9^( m - n)(e^ • e m _ n )cJ 1 Cn(di 1!/ - d mv ). 
The transformed Hamiltonian is 

H = e- s He s = ^ owc^' + Yj \ d ™ d ™ + \) + ( 118 ) 

n^n' ma ^ ' 

^ V(j\ - n'^CnC^Cn' - cl.Cn- 
n^n' n 

The last term describes the energy gained by polarons due to e-ph interaction. 
E p is the familiar polaron level shift, 



E p =w^^(m-n)(e l/ -e m _ n ) 2 , 



(119) 
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which is independent of n. The third term on the right-hand side in Eq.(130) 
is the polaron-polaron interaction: 



where 



w(n - n') = V c (n - n') - V ph (n - n'), (120) 
V ph (n-ri) = 211)0^2 g v (m-n)g v (m-ri) x 

The phonon-induced interaction V p h is due to displacements of common ions 
by two electrons. The transformed hopping operator <r nn ' in the first term in 
Eq.(118) is given by 



ow = T(n - n') exp 



^2 [9A m - n )( e » ■ e m-n) 



(121) 



g u (m - n')(e„ • e m _ n ')] (<4, Q - d ma )] ■ 



This term perturbation at large A. Here we consider a particular lattice struc- 
ture (ladder), where bipolarons tunnell already in the first order in T(n), so 
that cr nn ' can be averaged over phonons. When T<wq the result is 



i(n - n') = ((a nn ,)) ph = T(n - n') exp[-. 9 2 (n - n')], 



(122) 



g 2 (n- n') = y]ff„(m- n)(e„ • e m _ n ) x 

[g v {m - n)(e„ • e m _ n ) - g v {m - n')(e„ ■ e m _ n /)] . 
The mass renormalization exponent can be expressed via E p and V p h as 



ff 2 (n-n') 
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UJQ 



E p -±V ph (n-n>) 



(123) 



Now phonons are "integrated out" and the polaronic Hamiltonian is given by 

H p = Ho + H pert , (124) 

H = -E p ^ 4 c n + ^ V(ll - Il'yicnC^Cn' , 

ii n^n' 

Hpert = ^ t(n- n')clcn>. 
njta' 

When V p h exceeds V c the full interaction becomes negative and polarons form 
pairs. The real space representation allows us to elaborate more physics behind 
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Fig. 4. Mechanism of the polaron-polaron interaction, (a) Together, two polarons 
(solid circles) deform the lattice more effectively than separately. An effective at- 
traction occurs when the angle <j) between xi and X2 is less than n/2. (b) A mixed 
situation: ion 1 results in repulsion between two polarons while ion 2 results in 
attraction. 



the lattice sums in Eq.(119) and Eq.(120) [80]. When a carrier (electron or 
hole) acts on an ion with a force f , it displaces the ion by some vector x = 
f/fc. Here k is the ion's force constant. The total energy of the carrier-ion 
pair is — f 2 /(2fc). This is precisely the summand in Eq.(119) expressed via 
dimcnsionless coupling constants. Now consider two carriers interacting with 
the same ion, Fig. 4a. The ion displacement is x = (fi + f2)/fc and the energy 
is — f^/(2fc) — f|/(2fc) — (fi ■ f2)/fc. Here the last term should be interpreted as 
an ion-mediated interaction between the two carriers. It depends on the scalar 
product of fi and f 2 and consequently on the relative positions of the carriers 
with respect to the ion. If the ion is an isotropic harmonic oscillator, as we 
assume here, then the following simple rule applies. If the angle <f> between fi 
and f 2 is less than 7r/2 the polaron-polaron interaction will be attractive, if 
otherwise it will be repulsive. In general, some ions will generate attraction, 
and some repulsion between polarons, Fig. 4b. 

The overall sign and magnitude of the interaction is given by the lattice 
sum in Eq.(120), the evaluation of which is elementary. One should also note 
that according to Eq.(123) an attractive e-ph interaction reduces the polaron 
mass (and consequently the bipolaron mass) , while repulsive e-ph interaction 
enhances the mass. Thus, the long-range nature of the e-ph interaction serves 
a double purpose. Firstly, it generates an additional inter-polaron attraction 
because the distant ions have small angle <f>. This additional attraction helps 
to overcome the direct Coulomb repulsion between polarons. And secondly, 
the Frohlich interaction makes the bipolarons lighter. 
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The many-particle ground state of H depends on the sign of the polaron- 
polaron interaction, the carrier density, and the lattice geometry. Here we 
consider the zig-zag ladder, Fig. 5a, assuming that all sites are isotropic two- 
dimensional harmonic oscillators. For simplicity, we also adopt the nearest- 
neighbour approximation for both interactions, g v i\) = g, V c (n) = V c , and for 
the hopping integrals, T(m) = Tmn for / = n = m = a, and zero otherwise. 
Hereafter we set the lattice period a = 1. There are four nearest neighbours 
in the ladder, z = 4. Then the one-particle polaronic Hamiltonian takes the 
form 

H P = ~ E p X)(4cn + pIpu) + (125) 

n 

XX( C ™+1 C ™ +Pi+iP™) + t(pt c n + pl-lCn) + H.C.], 
n 

where c„ and p n are polaron annihilation operators on the lower and upper 
legs of the ladder, respectively, Fig. 5b. Using Eqs.(119), (120) and (123) we 
find 



E p = 4 3 2 w , (126) 

The Fourier transform of Eq.(125) into momentum space yields 

H p = ^2(2t' cos k-E p ){clc k +plp k )+ (127) 
fc 

t^ + e lk ) P {c k + H.c.}. 
fc 

A linear transformation of c k and p k diagonalises the Hamiltonian. There are 
two overlapping polaronic bands, 

E x (k) = -£ p + 2t'cos(fc)±2icos(fc/2) (128) 

with the effective mass m* — 2/\At' ± t\ near their edges. 

Let us now place two polarons on the ladder. The nearest neighbour inter- 
action is v = V c — Ep/2, if two polarons are on different legs of the ladder, and 
v = V c — Ep/4:, if both polarons are on the same leg. The attractive interaction 
is provided via the displacement of the lattice sites, which are the common 
nearest neighbours to both polarons. There are two such nearest neighbours 
for the intersite bipolaron of type A or B, Fig. 5c, but there is only one com- 
mon nearest neighbour for bipolaron C, Fig.5d. When V c > E p /2, there are 
no bound states and the multi-polaron system is a one-dimensional Luttingcr 
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Fig. 5. One-dimensional zig-zag ladder, (a) Initial ladder with the bare hopping 
amplitude T(a). (b) Two types of polarons with their respective deformations, (c) 
Two degenerate bipolaron configurations A and B. (d) A different bipolaron config- 
uration, C, whose energy is higher than that of A and B. 



liquid. However, when V c < E p /2 and consequently v < 0, the two polarons 
are bound into an inter-site bipolaron of types A or B. 

It is quite remarkable that the bipolaron tunnelling in the ladder already 
appears in the first order with respect to a single-electron tunnelling. This case 
is different from both onsite bipolarons discussed above, and from intersite 
chain bipolarons of Ref. [81] , where the bipolaron tunnelling was of the second 
order in T(a). Indeed, the lowest energy degenerate configurations A and 
B are degenerate. They are coupled by H pert . Neglecting all higher-energy 
configurations, we can project the Hamiltonian onto the subspace containing 
A, B, and empty sites. 

The result of such a projection is the bipolaronic Hamiltonian 

H b = (v c - \e v \ J2[A{A n +BlB n }-t'J2[BiA n +Bl_ 1 A n + H.c.}, (129) 

^ ^ n n 

where A n = c n p n and B n — p n c n+ \ are intersite bipolaron annihilation op- 
erators, and the bipolaron-bipolaron interaction is dropped (see below). The 
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Fourier transform of Eq.(129) yields two bipolaron bands, 

E 2 (k) = V c - -E p ± 2i'cos(fc/2). (130) 

with a combined width 4|t'|. The bipolaron binding energy in zero order with 
respect to t, t' is 

A = 2Bk (0) -E 2 (0) = ^--V c . (131) 
The bipolaron mass near the bottom of the lowest band, to** = 2/t', is 



to** = 4m* 



1 + 0.25exp ( 

8uj 



(132) 



The numerical coefficient 1/8 in the exponent ensures that to** remains of the 
order of to* even at large E p . This fact combines with a weaker renormalization 
of to* providing a superlight bipolaron. 

In models with strong intersite attraction there is a possibility of clusteri- 
sation. Similar to the two-particle case above, the lowest energy of n polarons 
placed on the nearest neighbours of the ladder is found as 

(\jT — 1 

E n = (2n - 3)V C —E p (133) 

for any n > 3. There are no resonating states for a n-polaron configuration 
if n > 3. Therefore there is no first-order kinetic energy contribution to their 
energy. E n should be compared with the energy Ei + (n — \)Eij1 of far 
separated (n — l)/2 bipolarons and a single polaron for odd n > 3, or with 
the energy of far separated n bipolarons for even n > 4. "Odd" clusters are 
stable if 

V c < E p , (134) 

6n-10 p v ' 

and "even" clusters are stable if 

n — 1 
on — 12 

As a result we find that bipolarons repel each other and single polarons at 
V c > 3E p /8. If Vc is less than 3E p /8 then immobile bound clusters of three and 
more polarons could form. One should notice that at distances much larger 
than the lattice constant the polaron-polaron interaction is always repulsive, 
and the formation of infinite clusters, stripes or strings is impossible (see 
also [47]). Combining the condition of bipolaron formation and that of the 
instability of larger clusters we obtain a window of parameters 

^E p < V c < \e p , (136) 

where the ladder is a bipolaronic conductor. Outside the window the ladder 
is either charge-segregated into finite-size clusters (small V c ), or it is a liquid 
of repulsive polarons (large V c ). 
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There is strong experimental evidence for superlight intersite bipolarons 
in cuprate superconductors (see below), where they form in-plane oxygen - 
apex oxygen pairs (so called apex bopolarons) and/or in- plane oxygen-oxygen 
pairs [301 H21 IBS] • While the long-range Frohlich interaction combined with 
Coulomb repulsion might cause clustering of polarons into finite-size quasi- 
metallic mesoscopic textures, the analytical [83 j and QMC [84] studies of 
mesoscopic textures with lattice deformations and Coulomb repulsion show 
that pairs dominate over phase separation since bipolarons effectively repel 
each other (see also [37] . 



8 Bipolaronic superconductivity 

In the subspace with no single polarons, the Hamiltonian of electrons strongly- 
coupled with phonons is reduced to the bipolaronic Hamiltonian written in 
terms of creation, b^ m = cj^ej^ and annihilation, b m , bipolaron operators as 

(137) 



i(m - m l )b\ a b ml + ^v{m - m')n m n 



where v(m — m') is the bipolaron-bipolaron interaction, n m — b^bm, and 
the position of the middle of the bipolaron band is taken as zero. There are 
additional spin quantum numbers S = 0, 1; S z = 0, ±1, which should be added 
to the definition of b m in the case of intersite bipolarons, which tunnel via the 
one-particle hopping. This "crab-like" tunnelling, Fig. 5, results in a bipolaron 
bandwidth of the same order as the polaron one. Keeping this in mind we can 
apply Hb, Eq.(137) to both on-site and/or inter-site bipolarons, and even to 
more extended non-overlapping pairs, implying that the site index m is the 
position of the centre of mass of a pair. 

Bipolarons are not perfect bosons. In the subspace of pairs and empty sites 
their operators commute as 

km&L + tf m b m = 1, (138) 
&m&L- 6 L' & m = ° ( 139 ) 

for m ^ m'. This makes useful the pseudospin analogy [4"T] . 

bl = S x m - iS« (140) 

and 

b f m b m = \-S z m (141) 

with the pseudospin 1/2 operators S x ' y ' z — \ti^,z- — 1/2 corresponds 
to an empty site m and = —1/2 to a site occupied by the bipolaron. 
Spin operators preserve the bosonic nature of bipolarons, when they are on 
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different sites, and their fermionic internal structure. Replacing bipolarons by 
spin operators we transform the bipolaronic Hamiltonian into the anisotropic 
Heisenberg Hamiltonian, 



(142) 



This Hamiltonian has been investigated in detail as a relevant form for mag- 
netism and also for quantum solids like a lattice model oi 4 He. However, while 
in those cases the magnetic field is an independent thermodynamic variable, 
in our case the total "magnetization" is fixed, 

^E« 5 m» = ^-^ (143) 
m 

if the bipolaron density rib is conserved. Spin 1/2 Heisenberg Hamiltonian, 
Eq.(142) cannot be solved analytically. Complicated commutation rules for 
bipolaron operators make the problem hard, but not in the limit of low atomic 
density of bipolarons, nj«l (for a complete phase diagram of bipolarons on a 
lattice see Refs. [41] [85]). In this limit we can reduce the problem to a charged 
Bose gas on a lattice [8 6) . Let us transform the bipolaronic Hamiltonian to a 
representation containing only the Bose operators a m and a^ m defined as 



V = ][>(aL) fe am + \ (144) 

k—0 

oo 

bl = Y / Mal) k+1 a k m , (145) 

k=0 

where 

a m o} m , - a? ml a m = 5 m . m >. (146) 

The first few coefficients (3k are found by substituting Eqs.(144) and (145) 
into Eqs.(138) and (139), 

1 /3 

/3 = l,/3i = -l,/3 2 = - + ^-. (147) 
We also introduce bipolaron and boson tf'-operators as 

«P(r) = ^=^(5(r-m)& m , (148) 



^ r ) = -^5>( r - m ) a ™- (149) 



The transformation of the field operators takes the form 
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<£(r) 



^(r^fr) (1/2 + 73/6)^(^^)^(1)^^) 



A 2 

Then we write the bipolaronic Hamiltonian as 



<P(r). 

(150) 



ff 6 = /" dr J dr'^(r)i(r - r')^(r') + H d + H h + , (151) 



where 



Hd = \ j dr j dr '^ T ~ r ')^ (r)^ (r'Pir'^ir), (152) 
is the dynamic part, 

H k = J dr J dr't(r - r') x (153) 
[i?t( r )#t( r ')^( r ')^( r ') + &\ T )ipi(r)W(T)&(r r )] . 

is the kinematic (hard-core) part due to the "imperfect" commutation rules, 
and ff( 3 ) includes three- and higher-body collisions. Here 



** ik-(r— r ) 



v(r-r') = ^Y,^^' 



(r-r') 



A 

k 

Ok = X)m^o ^( m ) exp(«k • m) is the Fourier component of the dynamic inter- 
action and 

e** = t(m) exp(-ik ■ m) (154) 

is the bipolaron band dispersion. contains powers of the field operator 
higher than four. In the dilute limit, rib ^ 1) only two-particle interactions are 
essential which include the short-range kinematic and direct density-density 
repulsions. Because v already has the short range part v^ 2 \ Eq.(113), the 
kinematic contribution can be included in the definition of v. As a result 
Hf> is reduced to the Hamiltonian of interacting hard-core charged bosons 
tunnelling in the narrow band. 

To describe electrodynamics of bipolarons we introduce the vector poten- 
tial A(r) using the so-called Peierls substitution [57], 

t(m - m') -> t(m - m ') e l2eA ( m )-( m - m ') j 

which is a fair approximation when the magnetic field is weak compared with 
the atomic field, eHa 2 << 1. It has the following form, 

t(r - r') t(r, r') = ]T etUeA^ ^ (155) 
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in real space. If the magnetic field is weak, we can expand e£* in the vicinity 
of k = to obtain 

t (r,rV-£±|^ 4 (r- r '), (156) 



where ^ 

is the inverse bipolaron mass. Here we assume a parabolic dispersion near the 
bottom of the band, e£* ~ k 2 , so that 

gtM -/^t (r) | [v±^! + At U (r)+ (158) 



dvdv'v(v - r')^\r)^ \r')V(r)V(r'), 



where we add the bipolaron chemical potential \i. We note that the bipolaron- 
bipolaron interaction is the Coulomb repulsion, v[v) ~ 1/(60?") at large dis- 
tances, and the hard-core repulsion is not important in the dilute limit. The 
Hamiltonian Eq.(159) describes the charged Bose gas with the effective boson 
mass m** and charge 2e, which is a superconductor [88] . 



9 Bipolaronic superconductivity in cuprates 

The fact that Helium-4 and its isotope Helium-3 are well known Bose and 
Fermi superfluids, respectively, with very different superfluid transition tem- 
peratures (T c = 2.YTK in A He and T c = 0.0026^ in 3 He) already kindles the 
view that high-temperature superconductivity might derive from preformed 
real-space charged bosons rather than in the BCS state with Cooper pairs, 
which are correlations in momentum space. The possibility of real-space pair- 
ing of carriers in cuprates, as opposed to Cooper pairing, has been the subject 
of much discussion. Some authors dismissed any real-space pairs even in un- 
derdoped cuprates, where a low density of carriers appears to favor individual 
pairing rather than Cooper pairing. But on the other hand real-space pairing 
is strongly supported by our strong-coupling extension of the BCS theory since 
the e-ph interaction is very strong in cuprates. Also on the experimental side a 
growing number of other independent observations point to the possibility that 
high-T c superconductors may not be conventional Bardeen-Cooper-Schrieffer 
(BCS) superconductors, but rather derive from the Bose-Einstein condensa- 
tion (BEC) of real-space superlight small bipolarons. There is strong evidence 
for real-space pairing and the three-dimensional BEC in cuprates from un- 
usual upper critical fields [5§] and isotope effects [106j predicted by us and 
the A-like electronic specific heat [90] , parameter- free fitting of experimental 
T c with BEC T c [91], normal state pseudogaps [92] [93] and anisotropy [94], and 
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more recently from normal state diamagnetism [95j , the Hall-Lorenz numbers 
[Ml |9"7] , the normal state Nernst effect [5H] |M] , and the giant proximity effect 
(GPE) [100] , Here I briefly discuss a few of these remarkable observations (for 
more details see Ref. [33] and Part IV). 

9.1 Upper critical field, the Hall-Lorenz number, and isotope 
effects 

Magnetotransport [101] and thermal magnetotransport ,102, 96] data strongly 
support preformed pairs in cuprates. In particular, many high magnetic field 
studies revealed a non-BCS upward curvature of the upper critical field, 
H<a(T) and its non- linear temperature dependence in the vicinity of T c in 
a number of cuprates as well as in a few other unconventional superconduc- 
tors, Fig. 6. If unconventional superconductors are in the "bosonic" limit of 
preformed real-space pairs, such unusual critical fields are expected in accor- 
dance with the theoretical prediction for the Bose-Einstein condensation of 
charged bosons in an external magnetic field [89 . 

Notwithstanding, some " direct" evidence for the existence of a charge 2e 
Bose liquid in the normal state of cuprates is highly desirable. Alexandrov and 
Mott [103] discussed the thermal conductivity k; the contribution from the 
carriers given by the Wiedemann- Franz ratio depends strongly on the elemen- 
tary charge as ~ (e*)~ 2 and should be significantly suppressed in the case of 
e* = 2e compared with the Fermi-liquid contribution. As a result, the Lorenz 
number, L — (e/kB) 2 n e /(Ta) differs significantly from the Sommerfeld value 
L e = 7r 2 /3 of the standard Fermi- liquid theory, if carriers are double-charged 
bosons. Here K e , <r, and e are the electronic thermal conductivity, the electri- 
cal conductivity, and the elementary charge, respectively. Ref. 103J predicted 
a rather low Lorenz number for bipolarons, L = 6L e / (Air 2 ) « 0.15L e , due to 
the double charge of carriers, and also due to their nearly classical distribution 
function above T c . 

The extraction of the electron thermal conductivity has proven difficult 
since both the electron term, n e and the phonon term, K p h are comparable to 
each other in the cuprates. A new way to determine the Lorenz number has 
been realized by Zhang et al. [102] . based on the thermal Hall conductivity. 
The thermal Hall effect allowed for an efficient way to separate the phonon 
heat current even when it is dominant. As a result, the "Hall" Lorenz number, 
Lh = {e/kB) K xy /(Ta xy ), has been directly measured in YBa,2Cu?,0§.^ be- 
cause transverse thermal n xy and electrical a xy conductivities involve only the 
electrons. Remarkably, the measured value of Lh just above T c is about the 
same as predicted by the bipolaron model, Lh ~ 0.15L e . The experimental 
Lh showed a strong temperature dependence, which violates the Wiedemann- 
Franz law. This experimental observation has been accounted for by taking 
into account thermally excited polarons and also triplet pairs in the bipolaron 
model [US], Fig.7. 
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Fig. 6. Resistive upper critical field [101] (determined at 50% of the tran- 
sition) of cuprates, spin-ladders and organic superconductors scaled accord- 
ing to the Bose-Einsten condensation field of charged bosons [89], H C 2(T) cc 
[fo(f — t)/t+ 1 — t 1 / 2 ] 3 2 with t = T/T c . The parameter b is proportional to the 
number of delocalised bosons at zero temperature, b is f (solid line), 0.02 (dashed- 
dotted line), 0.0012 (dotted line), and (dashed line). The inset shows a universal 
scaling of the same data near T c on the logarithmic scale. Symbols correspond to 
Ti-2201(»), LaimSro.isCuOiiA), £i-2201(x), Sz-2212(*), YBa 2 Cu 3 6+x {o), 
La2-xCe x CuC>4-y (squares), SY2Cai2Cu24 04i(+), and Bechgaard salt organic su- 
perconductor (V) . 

Another compelling evidence for (bi)polaronic carries in novel supercon- 
ductors was provided by the discovery of substantial isotope effects on T c and 
on the carrier mass [1041 1105] . The advances in the fabrication of the iso- 
tope substituted samples made it possible to measure a sizable isotope effect , 
a = —d\nT c /d\nM in many high-T c oxides. This led to a general conclusion 
that phonons are relevant for high T c . Moreover the isotope effect in cuprates 
was found to be quite different from the BCS prediction, a = 0.5 (or less). 
Several compounds showed a > 0.5 , and sometimes negative values of a were 
observed. 
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Fig. 7. The Hall Lorenz number Lh [SS] of charged bosons fits the experiment in 
YBa2Cu30e.95 102 . The pseudogap is taken as 675 K. The inset gives the ratio of 
the Hall Lorenz number to the Lorenz number in the model. 

Essential features of the isotope effect, in particular large values in low T c 
cuprates, an overal trend to lower value as T c increases, and a small or even 
negative a in some high T c cuprates can be understood in the framework of the 
bipolaron theory [106j . With increasing ion mass the bipolaron mass increases 
and the Bose-Einstein condensation temperature T c oc 1/m** decreases in the 
bipolaronic superconductor (section 8). On the contrary in polaronic super- 
conductors (section 6) an increase of the ion mass leads to a band narrowing 
enhancing the polaron density of states and increasing T c . Hence the isotope 
exponent of T c can distinguish the BCS like polaronic superconductivity with 
a < , and the Bose-Einstein condensation of small bipolarons with a > 0. 
Moreover, underdoped cuprates, which are certainly in the BEC regime, could 
have a > 0.5, as observed. 

The isotope effect on T c is linked with the isotope effect on the carrier 
mass, a m *, as |106] 

a = -d\nT c /d\nM = a m , [1 - Z/(X - (i c )], (159) 

where a m * = <i In to* /<i In M and Z = m/m* <C 1. In ordinary metals, where 
the Migdal approximation is believed to be valid, the renormalized effective 
mass of electrons is independent of the ion mass M because the electron- 
phonon interaction constant A does not depend on M. However, when the 
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e-ph interaction is sufficiently strong, the electrons form polarons dressed by 
lattice distortions, with an effective mass m* = mexp^Ep/oj). While E p in 
the above expression does not depend on the ion mass, the phonon frequency 
does. As a result, there is a large isotope effect on the carrier mass in pola- 
ronic conductors, a m * = (1/2) ln(m*/m) [106j . in contrast to the zero isotope 
effect in ordinary metals. Such an effect was observed in cuprates in the Lon- 
don penetration depth Xh of isotope-substituted samples [104] . The carrier 
density is unchanged with the isotope substitution of O 16 by O 18 , so that the 
isotope effect on \h measures directly the isotope effect on the carrier mass. 
In particular, the carrier mass isotope exponent a m * was found as large as 
a m » = 0.8 in Lai.gQsSro.wsCuOj,. 

More recent high resolution angle resolved photocmission spectroscopy 
|107j provided another compelling evidence for a strong e-ph interaction in 
cuprates. It revealed a fine phonon structure in the electron self-energy of 
the underdoped La2- x Sr x Cu04 samples. Remarkably, an isotope effect on 
the electron spectral function in Bi-2212 [108] has been discovered. These 
experiments together with a number of earlier optical |1091 11101 lllll 11121 11131 
1114) and neutron-scattering [115] experimental and theoretical studies firmly 
established the strong coupling of carries with optical phonons in cuprates 
(see also Part IV). 

9.2 Normal state diamagnetism: BEC versus phase fluctuations 

Above T c the charged bipolaronic Bose liquid is non-degenerate and below T c 
phase coherence (ODLRO) of the preformed bosons sets in. The state above 
T c is perfectly "normal" in the sense that the off-diagonal order parameter 
(i.e. the Bogoliubov-Gor'kov anomalous average jF(r, r') = (^(r)-0f (r')) is 
zero above the resistive transition temperature T c as in the BCS theory. Here 
ipl^(v) annihilates electrons with spin J,, f at point r. 

However in contrast with the bipolaron and BCS theories a significant 
fraction of research in the field of cuprate superconductors suggests a so-called 
phase fluctuation scenario [1161 11171 I118j . where !F(r, r') remains nonzero 
well above T c . I believe that the phase fluctuation scenario is impossible to 
reconcile with the extremely sharp resistive transitions at T c in high-quality 
underdoped, optimally doped and overdoped cuprates. For example, the in- 
plane and out-of-plane resistivity of Bi-2212, where the anomalous Nernst 
signal has been measured |117j . is perfectly "normal" above T c , Fig. 8, showing 
only a few percent positive or negative magnetoresistance [119] , explained with 
bipolarons [120] . 

Both in-plane [mj [1221 H23 H2H [125] and out-of-plane [Ell EH [T28] 
resistive transitions of high-quality samples remain sharp in the magnetic 
field providing a reliable determination of the genuine H C 2 (T) . The preformed 
Cooper-pair (or phase fluctuation) model [116] is incompatible with a great 
number of thermodynamic, magnetic, and kinetic measurements, which show 
that only holes (density x) 1 doped into a parent insulator are carriers both 
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in the normal and the superconducting states of cuprates. The assumption 
|116j that the superfluid density x is small compared with the normal-state 
carrier density is also inconsistent with the theorem [129] . which proves that 
the number of supercarriers at T = OK should be the same as the number of 
normal-state carriers in any clean superfluid. 

The normal state diamagnetism of cuprates provides another clear evi- 
dence for BEC rather than for the phase fluctuation scenario. A number of 
experiments (see, for example, [1301 11311 11321 11331 1 1341 1118] and references 
therein), including torque magnetometries, showed enhanced diamagnetism 
above T Cl which has been explained as the fluctuation diamagnetism in quasi- 
2D superconducting cuprates (see, for example Ref. 132J ) . The data taken 
at relatively low magnetic fields (typically below 5 Tcsla) revealed a cross- 
ing point in the magnetization M{T, B) of most anisotropic cuprates (e.g. 
Bi - 2212), or in M (T, B)/B 1 / 2 of less anisotropic YBCO [HE]. The depen- 
dence of magnetization (or M/B 1 / 2 ) on the magnetic field has been shown to 
vanish at some characteristic temperature below T c . However the data taken 
in high magnetic fields (up to 30 Tesla) have shown that the crossing point, 
anticipated for low-dimensional superconductors and associated with super- 
conducting fluctuations, does not explicitly exist in magnetic fields above 5 
Tesla [133] . 

Most surprisingly the torque magnetometery [130| 1133] uncovered a dia- 
magnetic signal somewhat above T c which increases in magnitude with applied 
magnetic field. It has been linked with the Nernst signal and mobile vortexes 
in the normal state of cuprates 118 . However, apart from the inconsistences 
mentioned above, the vortex scenario of the normal-state diamagnetism is 
internally inconsistent. Accepting the vortex scenario and fitting the magne- 
tization data in Bi — 2212 with the conventional logarithmic field dependence 
|118] , one obtains surprisingly high upper critical fields H c2 > 120 Tesla and a 
very large Ginzburg-Landau parameter, n — A/£ > 450 even at temperatures 
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Fig. 9. Diamagnetism of optimally doped Bi-2212 (symbols) 118 compared with 
magnetization of CBG [55] near and above T c (lines). 



close to T c . The in-plane low-temperature magnetic field penetration depth is 
A = 200 nm in optimally doped Bi — 2212 (see, for example |135j ). Hence the 
zero temperature coherence length £ turns out to be about the lattice con- 
stant, £ = 0.45nm, or even smaller. Such a small coherence length rules out the 
"preformed Cooper pairs" [116j . since the pairs are virtually not overlapped 
at any size of the Fermi surface in Bi — 2212 . Moreover the magnetic field 
dependence of M(T, B) at and above T c is entirely inconsistent with what 
one expects from a vortex liquid. While —M(B) decreases logarithmically at 
temperatures well below T c , the experimental curves [1301 11331 1118] clearly 
show that —M(B) increases with the field at and above T c , just opposite to 
what one could expect in the vortex liquid. This significant departure from 
the London liquid behavior clearly indicates that the vortex liquid does not 
appear above the resistive phase transition |130j . 




Some time ago we explained the anomalous diamagnetism in cuprates 
as the Landau normal-state diamagnetism of preformed bosons [136] . More 
recently the model has been extended to high magnetic fields taking into 
account the magnetic pair-breaking of singlet bipolarons and the anisotropy 
of the energy spectrum [95] . When the magnetic field is applied perpendic- 
ular to the copper-oxygen plains the quasi-2D bipolaron energy spectrum 
is quantized as E a = uj(n + 1/2) + 2£ c [l — cos(K z d)), where a comprises 
n = 0,1,2,... and in-plane K x and out-of-plane K z center-of-mass quasi- 
momenta, ui = 2eB / ^/m**m**, t c and d are the hopping integral and the 
lattice period perpendicular to the planes. We assume here that the spectrum 
consists of two degenerate branches, so-called "x" and "y" bipolarons as in 
the case of apex intersite pairs [40) with anisotropic in-plane bipolaron masses 
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m** = m and m** ~ 4m. Expanding the Bose- Einstein distribution function 
in powers of exp[(/i — E)/T] with the negative chemical potential /i one can 
after summation over n readily obtain the boson density 

n b = —r>Io(2t c rT) — , 160) 

7ra ' 1 — exp(— ujr/T) 

r—l x 11 

and the magnetization, 

M(T, B) = -n m + ^f> (^) x (161) 
exp[(/x-w/2-2t c )r/T] /l wexp(-wr/T) 



1 — cxp(— Lor/T) \r fcflT[l — exp(— ur/T)] 

Here /ib = e/ y/m**m** and Jo (^) is the modified Bessel function. At low tem- 
peratures T — > Schafroth's result [55J is recovered, M(0,B) = —rib^b- The 
magnetization of charged bosons is field-independent at low temperatures. At 
high temperatures, T>T C the chemical potential has a large magnitude, and 
we can keep only the terms with r = 1 in Eqs.(160,161) to obtain M(T, B) = 
—nb/ibtu I (6T) at T 3> T c 3> u>, which is the familiar Landau orbital diamag- 
netism of nondegenerate carriers. Here T c is the Bose-Einstein condensation 
temperature T c = 3.31(n b /2) 2 / 3 / {m* x *m* y *m* c *) 1 /' i , with m c = l/2\t c \d 2 . 

Comparing with experimental data one has to take into account a tem- 
perature and field depletion of singlets due to their thermal excitations 
into spin-split triplet states, nb(T,B) = n c [l — olt — (B/B*) 2 ]. Here a = 
3(2M)- 1 [^(e J/Tc -l)- 1 -T c ln(l--e- J / T 0], HbB* = (2T c n c t) 1 / 2 sinh( J/2T C ), 
[iB ~ 0.93 x 10 -23 Am 2 is the Bohr magneton, n c is the density of singlets 
at T = T c in zero field, r = T/T c — 1, J is the singlet-triplet exchange en- 
ergy, and 2t is the triplet bandwidth. As a result, Eq.(161) fits remarkably 
well the experimental curves in the critical region of optimally doped Bi-2212, 
Fig.9, with n c n b = 2100A/m, T c = 90K, a = 0.62 and B* = 56 Tesla, which 
corresponds to the singlet-triplet exchange energy J w 20K. 

On the other hand the experimental data, Fig.9, contradict BCS and the 
phase-fluctuation scenarios [1161 1118] . Indeed, if we define a critical expo- 
nent as S = hxB/ln\M(T,B)\ for B -> 0, the T dependence of 5(T) in 
the charged Bose gas (CBG) is dramatically different from the Berezinski- 
Kosterlitz-Thouless (BKT) transition critical exponents (as proposed in the 
phase fluctuation scenario), but it is very close to the experimental [118] S(T), 
Fig. 10. 

Also the large Nernst signal, allegedly supporting vortex liquid in the 
normal state of cuprates |118j . has been explained as the normal state phe- 
nomenon owing to a partial localization of charge carriers in a random po- 
tential inevitable in cuprates [55] . The coexistence of the large Nernst signal 
and the insulating-like resistivity in slightly doped cuprates sharply disagrees 
with the vortex scenario, but agrees remarkably well with our theory [99j . 
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9.3 Giant proximity effect 

Several groups reported that in the Josephson cuprate SNS junctions super- 
current can run through normal iV-barriers as thick as 100 nm in a strong 
conflict with the standard theoretical picture, if the barrier is made from 
non-superconducting cuprates. Using an advanced molecular beam epitaxy, 
Bozovic et al. |137j proved that this giant proximity effect (GPE) is intrin- 
sic, rather than extrinsic caused by any inhomogeneity of the barrier. Hence 
GPE defies the conventional explanation, which predicts that the critical cur- 
rent should exponentially decay with the characteristic length of about the 
coherence length, which is £ < 1 nm in the cuprates. 

This effect can be broadly understood as the Bose-Einstein condensate 
tunnelling into a cuprate semiconductor [100J. Indeed the chemical potential 
/i remains in the charge-transfer gap of doped cuprates like La2-xSr x Cu04 
|138j because of the bipolaron formation. The condensate wave function, ip(Z), 
is described by the Gross-Pitaevskii (GP) equation. In the superconducting 
region, Z < 0, near the SN boundary, Fig. 11, the equation is 

' <,lr[Z] -[V\^Z)?-mZ\ (162) 



2m** dZ 2 



where V is a short-range repulsion of bosons, and m** is the boson mass along 
Z. Deep inside the superconductor \i[)(Z)\ 2 = n s and /i = Vn s , where the 
condensate density n s is about x/2, if the temperature is well below T c of the 
superconducting electrode (the in-plane lattice constant a and the unit cell 
volume are taken as unity). 
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The normal barrier at Z > is an underdoped cuprate semiconductor 
above its transition temperature, where the chemical potential fj, lies below 
the bosonic band by some energy e, Fig. IT. For quasi-two dimensional bosons 
one readily obtains 

e{T) = -rin(l - e- To/T ), (163) 

where Tq — nx /to** , m** is the in-plane boson mass, and x' < x is the doping 
level of the barrier. Then the GP equation in the barrier can be written as 



d 2 ^(Z) 



2m** dZ 2 



= [V\^(Z)\ 2 + e}^(Z). (164) 



Introducing the bulk coherence length, £ = l/(2m**n s V) 1 ^ 2 and dimension- 

1/2 

less f(z) = ip(Z)/n s , p. — e/n s V, and z — Z/£, one obtains for a real 
/(*) 

S = f - /. (165) 



if z < 0, and 

/ 3 + /i/, (166) 



dz 2 



if z > 0. These equations can be readily solved using first integrals of motion 
respecting the boundary conditions, /(— oo) = 1, and /(oo) = 0, 

| = -(l/2 + / 4 /2-/ 2 ) 1 /2, (i 67 ) 

and 

f = -(/V 2 + / 4 /2) 1/2 , (168) 
dz 

for 2 < and z > 0, respectively. The solution in the superconducting elec- 
trode is given by 



f(z) = tanh 



-2-.'* = + 0.51„ 2 " 2 < 1+ <'> 1/2 + 1 



2 1 /2(l+/i)l/2_l 



(169) 



It decays in the close vicinity of the barrier from 1 to /(0) = [2(1 + /i)] -1 / 2 in 
the interval about the coherence length £. On the other side of the boundary, 
z > 0, it is given by 

f( z ) = ( 2 A) 1/2 (17Q) 

n > sinh{ 2 /i 1 /2 + i n [2( / i(l + ^))i/2 + (i + 4 / i(l + /i))i/2]}- ^ > 

Its profile is shown in Fig. 11. Remarkably, the order parameter penetrates 
into the normal layer up to the length Z* ss (A) _1 ^ 2 £, which could be larger 
than £ by many orders of magnitude, if jl is small. It is indeed the case, if 
the barrier layer is sufficiently doped. For example, taking x' = 0.1, c-axis 
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Fig. 11. BEC order parameter at the SN boundary for fl = 1.0,0.1,0.01 and 
< 0.001 (upper curve). 

m** = 2000m e , in-planc to** = 10to 6 [23] , a = 0.4 nm, and £ = 0.6 nm, yields 
T Q « 140 K and (A) 1/2 ~ 5000 at T = 10K. Hence the order parameter could 
penetrate into the normal cuprate semiconductor up to more than a thousand 
coherence lengths as observed |137j . If the thickness of the barrier L is small 
compared with Z*, and (/x) 1 / 2 <C 1, the order parameter decays following the 
power law, rather than exponentially, 



Hence, for L < Z* , the critical current should also decay following the power 
law [100] , On the other hand, for an undoped barrier jl becomes larger than 
unity, fl oc ln(TO**T/7rx') — > oo for any finite temperature T when x 1 — * 0, 
and the current should exponentially decay with the characteristic length 
smaller that £, as is experimentally observed as well j!38j . As a result the 
bipolaron theory accounts for the giant and nil proximity effects in slightly 
doped semiconducting and undoped insulating cuprates, respectively. It pre- 
dicts the occurrence of a new length scale, %/ ^/2m**e(T), and explains the 
temperature dependence of the critical current of SNS junctions [100] , 

10 Conclusion 



s/2 



(171) 



z + 2 



Extending the BCS theory towards the strong interaction between electrons 
and ion vibrations, a charged Bose gas of tightly bound small bipolarons was 
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predicted by us [41j with a further prediction that high T c should exist in 
the crossover region of the e-ph interaction strength from the BCS-like to 
bipolaronic superconductivity [18] . 

For very strong electron-phonon coupling, polarons become self-trapped 
on a single lattice site. The energy of the resulting small polaron is given as 
— E p = —Xzt. Expanding about the atomic limit in hopping integrals t (which 
is small compared to E p in the small polaron regime, A > 1) the polaron mass 
is computed as m* — mo exp^Ep/wo) , where loq is the frequency of Einstein 
phonons, mo is the rigid band mass on a cubic lattice, and 7 is a numerical 
constant. For the Holstcin model, which is purely site local, 7 = 1. Bipolarons 
are on-site singlets in the Holstein model and their mass m*^ appears only 
in the second order of t [41] scaling as m*^ oc (m*) 2 for uj ^> A , and as 
m*H oc (m*) 4 in a more realistic regime u> <C A (section 7). Here A = 2E p — U 
is the bipolaron binding energy, and U is the on-site (Hubbard) repulsion. 
Since the Hubbard U is about 1 eV or larger in strongly correlated materials, 
the electron-phonon coupling must be large to stabilize on-site bipolarons and 
the Holstein bipolaron mass appears very large, mjf/mo > 1000, for realistic 
values of the phonon frequency. 

This estimate led some authors to the conclusion that the formation of 
itinerant small polarons and bipolarons in real materials is unlikely [140j, and 
high-temperature bipolaronic superconductivity is impossible |141j . However, 
one should note that the Holstein model is an extreme polaron model, and 
typically yields the highest possible value of the (bi)polaron mass in the strong 
coupling limit. Many advanced materials with low density of free carriers and 
poor mobility (at least in one direction) are characterized by poor screening 
of high-frequency optical phonons and are more appropriately described by 
the long-range Frohlich electron-phonon interaction [40]. For this interaction 
the parameter 7 is less than 1 (7 pa 0.3 on the square lattice and 7 p=> 0.2 on 
the triangular lattice), reflecting the fact that in a hopping event the lattice 
deformation is partially pre-existent. Hence the unscreened Frohlich electron- 
phonon interaction provides relatively light small polarons, which are several 
orders of magnitude lighter than small Holstcin polarons. 

As shown above FCM is reduced to an extended Hubbard model with inter- 
site attraction and suppressed double-occupancy in the limit of high phonon 
frequency u> > t and large on-site Coulomb repulsion. Then the Hamiltonian 
can be projected onto the subspace of nearest neighbor intersite bipolarons. 
In contrast with the crawler motion of on-site bipolaron, the intersite bipo- 
laron tunnelling is a crab- like, so that its mass scales linearly with the polaron 
mass (m** p=> 4m* on the staggered chain [50]) as confirmed numerically us- 
ing CTQMC algorithm by Kornilovitch [32] • As a result , the crab bipolarons 
could bose-condense already at the room temperature [70], 

We believe that the following recipe is worth investigating to look for 
room-temperature superconductivity [70] : (a) The parent compound should 
be an ionic insulator with light ions to form high-frequency optical phonons, 
(b) The structure should be quasi two-dimensional to ensure poor screening of 
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high-frequency c-axis polarized phonons, (c) A triangular lattice is desirable 
in combination with strong, on-site Coulomb repulsion to form the superlight 
crab bipolaron, and (d) Moderate carrier densities are required to keep the sys- 
tem of small bipolarons close to the dilute regime. I believe that most of these 
conditions are already met in cuprate superconductors. As discussed above 
there is strong evidence for 3D bipolaronic BEC in cuprates from unusual 
upper critical fields and the electronic specific heat, normal state pseudogaps 
and anisotropy, normal state diamagnetism, the Hall-Lorenz numbers, and 
the giant proximity effect. 
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